#include "DDVAnalysis.hh"
#include "DDVData.hh"
#include "DDChargeAnalysis.hh"
#include "DDEventManager.hh"
#include "TMath.h"
#include <iostream>
#include <set>
#include <cmath>

using namespace std;

int DDChargeAnalysis::Analyze(DDVData* data)
{
	int event;
	float min_num_electrons;
	float max_num_electrons;
	float total_num_electrons;

	// Loop over channels
	//map<int, TH1F*>::iterator channel;
	set<int>::iterator channel;
	//for (channel=data->GetCurrentTrace().begin();channel!=data->GetCurrentTrace().end(); ++channel)
	for (channel=(*(data->get_header()->get_channels())).begin();channel!=(*(data->get_header()->get_channels())).end(); ++channel)
	{
		// Rangefinding loop
		for(event=0;event<data->get_header()->get_num_traces();event++)
		{
			cout << data->get_header()->get_num_traces() << endl;
			data->GoToEvent(event);
			minmax(data, *channel, event, &min_num_electrons, &max_num_electrons, &total_num_electrons); 
		}
	
		// Create histogram
		create_histo(data, *channel, &min_num_electrons, &max_num_electrons, &total_num_electrons);
	
		//Analysis loop
		for(event=0;event<data->get_header()->get_num_traces();event++)
		{
			data->GoToEvent(event);
			charge_histo(data, *channel); 
		}
	cout << "loopdone" << endl;
	}
	cout << "end loop" <<endl;
	return DDEventManager::kOK;
}

void DDChargeAnalysis::minmax(DDVData* data, int channel, int event, float* min_num_electrons, float* max_num_electrons, float* total_num_electrons)
{
	TH1F* trace = (data->GetCurrentTrace())[channel];

	cout << "mm pre int" << endl;

	float charge = float(trace->Integral())/50.0*data->get_header()->get_time_per_sample(channel);
	cout << "mm post int" << endl;

	cout << charge << endl;

	float num_electrons= charge/1.6e-19; //1.6e-19 is the charge of an electron

	if (event == 0)
	{
    	*max_num_electrons = num_electrons;
    	*min_num_electrons = num_electrons;
	}

    if (num_electrons > *max_num_electrons)
    {
    	*max_num_electrons = num_electrons;
    }

    if (num_electrons <  *min_num_electrons)
    {
    	*min_num_electrons = num_electrons;
    }

	//total_num_electrons+=num_electrons; 
}

void DDChargeAnalysis::create_histo(DDVData* data, int channel, float* min_num_electrons, float* max_num_electrons, float* total_num_electrons)
{
//	float maxrange = pow(10,ceil(log10(fabs(max_num_electrons)))) * (max_num_electrons/fabs(max_num_electrons));
//	float minrange = pow(10,floor(log10(fabs(min_num_electrons))))* (min_num_electrons/fabs(min_num_electrons));

	//float mean_num_electrons = total_num_electrons/num_traces_cut;
	float mean_num_electrons = *total_num_electrons/data->get_header()->get_num_traces();

	float maxrange = mean_num_electrons + 2*fabs(*max_num_electrons-mean_num_electrons);
	float minrange = mean_num_electrons - 2*fabs(*min_num_electrons-mean_num_electrons);

    current_histos[channel] = new TH1F("trace_electrons", "Electrons per trace", 100, minrange, maxrange);
}

void DDChargeAnalysis::charge_histo(DDVData* data, int channel){


  TH1F* trace = data->GetCurrentTrace()[channel];

  cout << endl;
  cout << "charge pre int" << endl;
  cout << endl;
  float charge = float(trace->Integral())/50.0*data->get_header()->get_time_per_sample(channel); //Coulombs, I=dq/dt V=(dq/dt)*R integral(V) dt / R = q, R from the oscilloscope = 50 ohms, division by timescale parameter to convert current from coulombs/timediv to coulombs/sec.
  cout << "charge post int" << endl;

  float num_electrons= charge/1.6e-19; //1.6e-19 is the charge of an electron

//  if (num_electrons < num_electrons_threshold)
//  {
//  	num_electrons_outside++;
//  }

	cout << "prefill" << endl;

  current_histos[channel]->Fill(num_electrons); // Fill the my_waveform electron histogram
	cout << "postfill" << endl;
  
//  if (writewaveform)
 //     my_waveform->Write(); // Write the waveform to the root file
	  //my_waveform->Smooth();
      //my_waveform->Write(); // Write the waveform to the root file

// ----Look at negative my_waveforms---
/*
      if (num_electrons < -5e27)
      {
//    writemy_waveform = true;
//    analysis->AnalyzeWaveform(my_data, time_per_div, &num_electrons, writemy_waveform);

      num_neg++;
//    my_waveform->SaveAs("
      }
*/
//-------------------------

}
/*
void Analysis::write_charge_histo()
{ 
	trace_electrons->Write();
	cout<< "---------------charge_histo-----------------"<<endl;
    cout<< max_num_electrons << " " << min_num_electrons << endl;
    cout<< "max num electrons" << max_num_electrons << endl;
    cout<< "min num electrons" << min_num_electrons << endl;
	cout<< "Num Neg: " << num_neg << endl;
	cout<< "Noise threshold (3 sigma below mean): " << num_electrons_threshold << endl;
	cout<< "Number of traces below threshold: " << num_electrons_outside << endl;
	cout<< "Percentage of traces below threshold: " << 100*(float(num_electrons_outside)/float(num_traces)) << " %" << endl;
	cout<< "------------------------------------------"<<endl;
	cout<< endl;
}
*/
